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ABSTRACT 

This paper concentrates on the flow prop- 
erties when one fluid displaces another fluid 
in a two-dimensional (2D) network of pores 
and throats. We consider the scale where 
individual pores enter the description and 
we use a network model to simulate the dis- 
placement process. 

We study the interplay between the pres- 
sure build up in the fluids and the dis- 
placement structure in drainage. We find 
that our network model properly describes 
the pressure buildup due to capillary and 
viscous forces and that there is good cor- 
respondence between the simulated evolu- 
tion of the fluid pressures and earlier results 
from experiments and simulations in slow 
drainage. 

We investigate the burst dynamics in 
drainage going from low to high injection 
rate at various fluid viscosities. The bursts 
are identified as pressure drops in the pres- 
sure signal across the system. We find that 
the statistical distribution of pressure drops 
scales according to other systems exhibiting 
self-organized criticality. We compare our 
results to corresponding experiments. 

We also study the stabilization mecha- 
nisms of the invasion front in horizontal 2D 
drainage. We focus on the process when 
the front stabilizes due to the viscous forces 
in the liquids. We find that the difference 
in capillary pressure between two different 
points along the front varies almost linearly 
as function of length separation in the di- 
rection of the displacement. The numer- 
ical results support new arguments about 
the displacement process from those earlier 
suggested for viscous stabilization. Our ar- 
guments are based on the observation that 



nonwetting fluid flows in loopless strands 
(paths) and we conclude that earlier sug- 
gested theories are not suitable to drainage 
when nonwetting strands dominate the dis- 
placement process. We also show that the 
arguments might influence the scaling be- 
havior of the front width as function of the 
injection rate and compare some of our re- 
sults to experimental work. 

1 INTRODUCTION 

Two-phase displacements in porous media 
have received much attention during the 
last two decades. In modern physics, the 
process is of great interest due to the variety 
of structures obtained when changing the 
fluid properties like wettability, interfacial 
tension, viscosities and displacement rate. 
The different structures obtained have been 
organized into three flow regimes: viscous 
fingering [1,2], stable displacement [3], and 
capillary fingering [4-6]. Viscous fingering 
is characterized by an unstable front of fin- 
gers that is generated when nonwetting and 
less viscous fluid is displacing wetting and 
more viscous fluid at relative high injection 
rate. The fingering structure is found to 
be fractal with fractal dimension D = 1.62 
[1,2]. Stable displacement is named after 
the relative flat and stable front that is be- 
ing generated when a nonwetting and more 
viscous fluid displaces a wetting and less 
viscous fluid at relative high injection rate. 
The last scenario, capillary fingering, is ob- 
tained when a nonwetting fluid very slowly 
displaces a wetting fluid. At sufficiently low 
injection rate the invasion fluid generates 
a pattern similar to the cluster formed by 
invasion percolation [4,7-9]. The displace- 
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merit is now solely controlled by the capil- 
lary pressure, that is the pressure difference 
between the two fluids across a meniscus in 
a pore. 

Fluid flow in porous media has also been 
intensively studied because of important 
applications in a wide range of different 
technologies. The most important areas 
that to a great extent depend on proper- 
ties of fluid flow in porous media, are oil 
recovery and hydrology. In oil recovery, 
petroleum engineers are continuously de- 
volving improved techniques to increase the 
amount of oil they are able to achieve from 
the oil reservoirs. In hydrology, one impor- 
tant concern is often to avoid pollution of 
ground water from human activity. 

The simulation model presented in this 
paper is developed to study the dynamics of 
the temporal evolution of the fluid pressures 
when a nonwetting fluid displaces a wetting 
fluid at constant injection rate. With the 
model we study the pressure in the fluids 
caused by the viscous forces as well as the 
capillary forces due to the menisci in the 
pores. The model porous medium consists 
of a tube network where the tubes are con- 
nected together to form a square lattice. 

Numerical simulations of fluid flow in 
porous media using a network of tubes 
was first proposed by Fatt [10] in 1956. 
Since then a large number of publications 
related to network models and pore-scale 
displacements have appeared in the litera- 
ture [1,3,11-23]. Often mentioned is the 
classic work of Lenormand et al. [3] who 
were the first to systematically classify the 
displacement structures into the three flow 
regimes: viscous fingering, stable displace- 
ment and capillary fingering. Including the 
work of Lenormand et al., it appears that 



most network models have been used to 
study statistical properties of the displace- 
ment structures or to calculate macroscopic 
properties like fluid saturations and relative 
permeabilities. 

There have been several attempts to sim- 
ulate the displacement process by using 
different types of growth algorithms. In 

1983 Wilkinson and Willemsen [9] formu- 
lated a new form of percolation theory, in- 
vasion percolation (IP), that corresponds to 
slow drainage, i.e. capillary fingering. In 

1984 Pater son [24] was the first to discover 
the remarkable parallels between diffusion- 
limited aggregation (DLA) [25] and viscous 
fingering. He also showed similarities be- 
tween anti-DLA and stable displacement. 
The disadvantage with the growth algo- 
rithms is that they do not contain any phys- 
ical time and they have so far not been suit- 
able to study the cross over between the 
different flow regimes. However, attempts 
have been made to use DLA and IP to study 
dynamics of viscous fingering [26] and slow 
drainage [27,28], respectively. 

In slow drainage it is observed that the 
invasion of nonwetting fluid occurs in a se- 
ries of bursts accompanied by sudden neg- 
ative drops in the pressure called Haines 
jumps [27-29] (see Fig. 1). This type of dy- 
namics is very important for the temporal 
evolution of the pressure in drainage, and in 
most network models the effect is neglected. 
Consequently, only very few network mod- 
els [21] have been used to study the inter- 
play between fluid pressures and displace- 
ment structures, and many questions ad- 
dressing this topic are still open. We will 
attempt to answer some of them in this 
paper, by making a model whose proper- 
ties are closer to those of real porous me- 
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Figure 1: Nonwetting fluid (white) invades 
a 2D porous medium initially filled with 
wetting fluid (shaded). As the nonwetting 
fluid is pumped into the system (a) the 
menisci move into narrower parts of the 
pore necks and the capillary pressure in- 
creases (b). During a burst the invading 
fluid covers new pores and the neighbor- 
ing menisci readjust back to larger radii and 
the capillary pressure decreases everywhere 
(c) [27]. The arrow in (a) is pointing at a 
pore neck having a shape of an hourglass. 



dia. To model the burst dynamics, we have 
been motivated by the hourglass shaped 
pore necks in Fig. 1. As a result we let 
the tubes in our network model behave as if 
they were hourglass shaped with respect to 
the capillary pressure. Thus, the capillary 
pressure of a meniscus starts at zero when 
the meniscus enters the tube and increases 
towards a maximum value at the middle of 
the tube where the tube is most narrow, be- 
fore the capillary pressure decreases to zero 
again when the meniscus leave the tube. 

The advantage of the above approach is 
a network model that reproduces the burst 
dynamics and the corresponding pressure 
evolution. We are also able to study in de- 
tail the capillary pressure of each meniscus 
along the front as it moves through the net- 
work. Similar measurements can hardly be 
done experimentally. 



We use the model to study the burst dy- 
namics going from low to high displacement 
rates. To do so, we examine the statistical 
properties of the sudden negative pressure 
drops due to bursts. We find that for a 
wide range of displacement rates and fluid 
viscosities, the pressure drops act in anal- 
ogy to theoretical predictions of systems 
exhibiting self-organized criticality, such as 
IP. Even at high injection rates, where the 
connection between the displacement pro- 
cess and IP is more open, the pressure drops 
behave similar to the case of extreme low 
injection rate, where IP apply. 

Further, we report on the behavior of 
the capillary pressure along the invasion 
front and investigate the stabilization mech- 
anisms of horizontal drainage. We present 
theoretical arguments predicting the behav- 
ior of the pressure along the front, and we 
conclude that the difference in pressure be- 
tween two different points along the front 
should depend almost linearly as function of 
the distance between the two points in the 
direction of the displacement. The theoret- 
ical arguments are based on the observation 
that the nonwetting fluid displaces the wet- 
ting fluid through separate loopless strands 
(paths). Numerical simulations of the cap- 
illary pressure along the front supports the 
theoretical arguments. We note that earlier 
suggested views [30-33] concerning the be- 
havior of the pressure along the front, is not 
compatible with our results. 

Unfortunately, the detailed modeling of 
the menisci's movements and their capillary 
pressures makes the model computationally 
heavy and reduces the system size that is 
attainable within feasible amount of CPU 
time. 

The paper is organized as follows. In Sec- 
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tion 2 we present the network model and 
describe briefly its numerical implementa- 
tion. The rest of the sections discuss the 
main results that we have got from numer- 
ical experiments with the network model. 
In section 3 we discuss the evolution of the 
simulated pressure in drainage and calcu- 
late the statistics of the bursts. We also 
compare some of our results to experimental 
work. Section 4 presents theoretical argu- 
ments about the stabilization mechanisms 
of the front in drainage. The arguments are 
supported by numerical simulations with 
the network model. This section also con- 
tains experimental work that is related to 
the numerical simulations. A summary of 
the results and concluding remarks are pro- 
vided in Section 5. 

2 SIMULATION 
MODEL 

The network model is thoroughly discussed 
in Ref. [34], and it has also been presented 
briefly in Refs. [35,36]. Therefore, only its 
main features are described in this section. 

The model porous medium consists of a 
square lattice of cylindrical tubes of length 
d oriented at 45° to the longest side of the 
lattice. Four tubes meet at each intersec- 
tion where we put a node having no volume. 
The disorder is introduced by (1) assigning 
the tubes a radius r chosen at random inside 
the interval [Xid, X 2 d] where < Ai < A 2 < 1 
or (2) moving the intersections a randomly 
chosen distance away from their initial po- 
sitions. The randomly chosen distances are 
less than 1/2 of the distance between the 
nearest neighbor intersections to avoid over- 




Figure 2: Example of a displacement struc- 
ture from one simulation. The nonwetting 
fluid (black) is injected from below and dis- 
places the wetting fluid (grey) that escapes 
along the top row. The front between the 
nonwetting and wetting phase is defined 
as the line separating the compact wetting 
phase and the nonwetting fluid. Note the 
trapped regions of wetting fluid that are left 
behind and surround by nonwetting fluid. 



lapping nodes. In (1) all tubes have equal 
d but different r. (2) results in a distorted 
square lattice giving the tubes different oTs. 
In (2) r = d/2a where a is the aspect ra- 
tio between the tube length and its radius. 
The reason for making a distorted lattice of 
tubes is to get closer to a real pore-throat 
geometry as shown in Fig. 1 [36]. 

Figure 2 shows an example of a displace- 
ment structure that is obtained from one 
simulation. The nonwetting fluid (black) 
of viscosity /i nw is injected along the in- 
let (bottom row) and displaces the wetting 
fluid (grey) of viscosity fi w . The fluids flow 
from the bottom to the top of the lattice, 
and there are periodical boundary condi- 
tions in the orthogonal direction. We as- 
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sume that the fluids are immiscible and in- 
compressible. 

A meniscus is located in the tubes where 
nonwetting and wetting fluids meet. The 
capillary pressure p c of a meniscus in a 
cylindrical tube of radius r is given by 
Young-Laplace law, 



P< 



2 ^ a 
= — cosy, 

r 



(1) 



under the assumption that the principal 
radii of the curvature of the meniscus are 
equal to the radius of the tube. 8 denotes 
the wetting angle between the cylinder wall 
and the wetting fluid, i.e. 0° < 9 < 90° in 
drainage. 

In the network model we treat the tubes 
as if they were hourglass shaped with re- 
spect to the capillary pressure. Therefore, 
we let the capillary pressure depend on 
where the meniscus is situated in the tube. 
Instead of Eq. (1) we let p c of a meniscus 
vary in the following way: 



27 

Pc = — 

r 



1 — cos (27r^| 



(2) 



Here we have assumed that the wetting fluid 
perfectly wets the medium, i.e. 9 = 0. In 
the above relation x denotes the position of 
the meniscus in the tube (0 < x < d), giv- 
ing that p c = at the entrance and at the 
exit of the tube and reaches a maximum of 
47/r in the middle of the tube (x = d/2). 
Practically, the wetting angle of a meniscus 
and thereby its capillary pressure may gen- 
erally be different depending on whether the 
meniscus retires from or invades the tube. 
To avoid numerical complications this effect 
is neglected in the present model. 

We solve the volume flux through each 
tube by using Hagen-Poiseuille flow for cy- 



lindrical tubes and the Washburn approxi- 
mation [37] for menisci under motion. Let 
qij denote the volume flux through the tube 
from the ith to the jth node, then we have 



Qij 



u 



P. 



(3) 



Here is the permeability of the tube 
(rfj/8) and is the cross section (nrfj) 
of the tube, /i^ denotes the effective vis- 
cosity, that is the sum of the volume frac- 
tions of each fluid inside the tube multiplied 
by their respective viscosities. The pres- 
sure drop across the tube is Apij = pj — pi, 
where pi and pj is the pressures at node i 
and j, respectively. The capillary pressure 
p c>i j is the sum of the capillary pressures 
of each menisci [given by Eq. (2)] that are 
present inside the tube. A tube partially 
filled with both liquids is allowed to con- 
tain at maximum two menisci. For a tube 
without menisci, p Ct ij = 0. We only con- 
sider horizontal flow, and therefore we ne- 
glect gravity. 

We have conservation of volume flux at 
each node giving 



2 Qa = °- 



(4) 



The summation on j runs over the near- 
est neighbor nodes to the ith node while i 
runs over all nodes that do not belong to 
the top or bottom rows, that is, the inter- 
nal nodes. Eqs. (3) and (4) constitute a set 
of linear equations which we solve for the 
nodal pressures p iy with the constraint that 
the pressures at the nodes belonging to the 
upper and lower rows are kept fixed. The 
set of equations is solved by using the Con- 
jugate Gradient method [38]. 
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In the simulations we impose the injec- 
tion rate Q on the inlet, therefore we have 
to find the pressure across the lattice AP, 
that corresponds to the given Q. Having 
found AP we use this pressure to calculate 
the correct piS in Eq. (3). In short, we find 
AP by considering the relation 



Q = AAP + B. 



(5) 



The first part of Eq. (5) results from Darcy's 
law for single phase flow through porous 
media. The second part comes from the 
capillary pressure between the two fluids 
(i.e. B = if no menisci are present in the 
network). Eq. (5) has two unknowns, A and 
B, which we calculate by solving Eq. (4) 
twice for two different applied pressures 
AP' and AP", across the lattice. From 
those two solutions we find the correspond- 
ing injection rates Q' and Q" . Inserting 
Q', Q", AP', and AP" into Eq. (5) results 
in two equations which we solve for A and 
B. Finally, we find the correct AP due to 
the imposed Q by rewriting Eq. (5), giving 
AP = (Q-B)/A. See Refs. [34,35] for fur- 
ther details on how the p^s are calculated 
after AP is found. 

Given the correct solution of pt we calcu- 
late the volume flux through each tube 
in the lattice, using Eq. (3). Having found 
the qi/s, we define a time step At such that 
every meniscus is allowed to travel at most 
a maximum step length Ax max during that 
time step. Each meniscus is moved a dis- 
tance (qij/(iij)At and the total time lapse 
is recorded before the nodal pressures pi, 
are solved for the new fluid configuration. 
Menisci that are moved out of a tube dur- 
ing a time step are spread into neighboring 
tubes as described in Refs. [34,35]. 



3 TEMPORAL EVO- 
LUTION OF FLUID 
PRESSURE 

To characterize the different fluid properties 
used in the simulations, we use the capil- 
lary number C a and the viscosity ratio M. 
The capillary number indicates the ratio be- 
tween viscous and capillary forces and in 
the simulations it is defined as 



a, = 



£ 7 



(6) 



Here Q is the injection rate of the nonwet- 
ting fluid, ji is the maximum viscosity of 
the nonwetting and wetting fluid and S is 
equal to the length of the inlet times the 
average thickness of the lattice, i.e. £ is the 
cross section of the inlet. 7 is the fluid- fluid 
interface tension. 

The viscosity ratio M, is defined as 



M = 



fir 



(7) 



where \i nw and [i w is the viscosity of the 
invading nonwetting fluid and the defending 
wetting fluid, respectively. 

3.1 TRAPPED FLUID AND 
PRESSURE BUILDUP 

The pressure across the system is found 
from Eq. (5) giving 



AP 



Q 

A 



P 



cgi 



where P cg = —B/A defines the global cap- 
illary pressure of the system. As will be- 
come clear below, P cg contains the capil- 
lary pressures of the menisci surrounding 
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Figure 3: AP (a), P cg (b), and A /A (c) as 
function of injection time. C a = 4.6 xl0~ 3 
and M = 100. The vertical dashed line is 
drawn at the saturation time, t s . 



the trapped wetting fluid (cluster menisci) 
and the capillary pressures of the menisci 
along the invasion front (front menisci) (see 
Fig- 2). 

Figure 3 shows the simulated pressures 
AP and P cg for a displacement at C a = 
4.6 x 10~ 3 and M = 100. The front width 
was observed to stabilize after some time t s , 
and a typical compact pattern of small clus- 
ters of wetting fluid developed behind the 
front. From Fig. 3 we observe that both AP 
and P cg increases as the more viscous fluid 
is pumped into the system. When t > t s 
they even tend to increase linearly as func- 
tion of time. 

The driving mechanism in the displace- 
ment is the pressure gradient between the 
inlet and the front causing a viscous drag 
on the trapped clusters. At moderate in- 
jection rates these clusters are immobile, 
thus the viscous drag is balanced by capil- 



lary forces along the interface of the cluster. 
On average the sum of the capillary forces 
from each cluster contributes to P cg by a 
certain amount making P cg proportional to 
the number of clusters behind the front. Af- 
ter the front has saturated with fully devel- 
oped clusters behind (t > t s ), the number 
of clusters are expected to increase linearly 
with the amount of injected fluid. Since 
the injection rate is held fixed we recognize 
that P cg must increase linearly as function 
of time. The argument does not apply when 
t < t s , due to the fractal development of the 
front before saturation. 

In Fig. 3 we have also plotted A /A which 
is the normalized difference between AP 
and P cg [see Eq. (8)]. A is equal to the 
proportionality factor between Q and AP 
when only one phase flows through the lat- 
tice (i.e. P cg = 0). We observe that A /A 
tends to increase linearly as function of time 
when t > t s . From Eq. (5) we interpret 
A as the total conductance of the lattice, 
and the reciprocal of that is the total re- 
sistance. The total resistance depends on 
the fluid configuration and the geometry 
of the network. Locally, the fluid config- 
uration changes as nonwetting fluid invades 
the system, however, the linear behavior 
of A /A indicates that the overall displace- 
ment structure is statistically invariant with 
respect to the injection time. That means, 
after the front has saturated (t > t s ) the 
displacement structure might be assigned a 
constant resistance per unit length. 

In the special case when M — 1 (vis- 
cosity matched fluids) the total resistance, 
1/A, was found to be constant indepen- 
dent of the injection rate or displacement 
structure. This somewhat surprising result 
might be explained by the following con- 
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Figure 4: P c f (a) and P cg (b) as function of 
injection time at C a = 3.5 xl0~ 4 and M = 
1.0 xlCT 3 . To avoid overlapping curves P cg 
was subtracted by 1000dyn/cm 2 before it 
was normalized. 



sideration. When M = 1 the effective vis- 
cosity fiij, of each tube is independent of 
the amount of wetting and nonwetting fluid 
that occupies the tube. Hence, each tube 
has a constant mobility of kij/fiij giving a 
constant total resistance of the network. 

At low C a we approach the regime of 
capillary fingering and the viscous drag on 
the clusters becomes negligible. Hence, P cg 
is no longer a linear function of the in- 
jection time, but reduces to that describ- 
ing the capillary pressure along the front. 
This is observed in Fig. 4 where we com- 
pare P cg with the calculated average cap- 
illary pressure along the front, P c f. In 
the simulations, P c f is calculated by taking 
the mean of the capillary pressures of the 
front menisci. From the figure we see that 
P cg ~ P c f, as expected. The big jumps in 
the pressure functions in Fig. 4 are caused 
by the variations in the capillary pressure 



as menisci move through the "hourglass" 
shaped tubes. The negative jumps are iden- 
tified as bursts where a meniscus proceeds 
abruptly [27, 29] to fill the tube with non- 
wetting fluid (see also Sec. 3.2 for further 
details). 

From the above discussion we conclude 
that the behavior of P cg at large times (t > 
t s ) may be formulated as 



Peg A mc H ~\~ Praf ■> 



(9) 



where A mc is the proportionality factor be- 
tween P cg and the average front position h 
above the inlet. A mc is given by the viscous 
drag on the clusters and P m f is the variation 
in the capillary pressure when the invasion 
front covers new tubes, h is only defined af- 
ter the front has saturated, i.e. h s < h < L, 
where h s is the average front position at t s 
and L is the length of the system. Since the 
injection rate is held fixed, h is proportional 
to the injection time t. In the limit of very 
low injection rates, A mc — > 0. 

When the average front position has 
reached the outlet, i.e. h = L in Eq. (9), 
only invading fluid flows through the sys- 
tem and P m f = 0. In this limit Darcy's law 
applied on the nonwetting phase gives U = 
(K e / fi nw )(AP/ L), where K e is the effec- 
tive permeability of the nonwetting phase. 
From Eqs. (8) and (9) we find that AP = 
Q/A + A mc L, which inserted into the Darcy 
equation gives 



l/a T + A mc /U' 



(10) 



Here a T = AL/Y, denotes the total con- 
ductivity of the lattice. Thus, we might 
consider the effective permeability of the 
nonwetting phase as a function of the con- 
ductivity of the lattice and an additional 
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term due to the viscous drag on the clus- 
ters (A mc /U). Note that the U dependency 
in Eq. (10) only indicates changes in A mc 
between displacements executed at different 
injection rates. The behavior when the flow 
rate changes during a given displacement is 
not discussed here. 

3.2 BURST DYNAMICS 

In the simulations a burst starts when the 
pressure drops suddenly and stops when 
the pressure has raised to a value above 
the pressure that initiated the burst (see 
Fig. 5). Thus, a burst may consist of a 
large pressure valley containing a hierarchi- 
cal structure of smaller pressure jumps {i.e. 
bursts) inside. A pressure jump, indicated 
as Ap in Fig. 5, is the pressure difference 
from the point when the pressure starts de- 
creasing minus the pressure when it stops 
decreasing. We define the size of the pres- 
sure valley (valley size) to be x = Hi Api, 
where the summation index % runs over all 
the pressure jumps Api inside the valley. 
The definition is motivated by experimen- 
tal work in Ref. [28]. For slow displace- 
ments we have that \ is proportional to the 
geometric burst size s, being invaded dur- 
ing the pressure valley. This statement has 
been justified in Ref. [28], where it was ob- 
served that in stable periods, the pressure 
increased linearly as function of the volume 
being injected into the system. Later, in an 
unstable period where the pressure drops 
abruptly due to a burst, this pressure drop 
is proportional to s. At fast displacements 
the pressure may no longer be a linear func- 
tion of the volume injected into the system. 
Therefore, a better estimate of s there, is 
to compute the time period T of the pres- 
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Figure 5: The pressure signal as function of 
injection time, P(t), for one simulation at 
low displacement rate in a narrow time in- 
terval. The horizontal line defines the pres- 
sure valley of a burst that last a time period 
T — t 2 — 1±. Note that the valley may con- 
tain a hierarchical structure of smaller val- 
leys inside. The vertical line indicates the 
size of a local pressure jump Ap inside the 
valley. 



sure valley (Fig. 5). Since the displacements 
are performed with constant rate, it is rea- 
sonable to assume that T is always propor- 
tional to the volume being injected during 
the valley and hence, T oc s. 

In Fig. 6 we have plotted the hierarchi- 
cal valley size distribution N a \\(x), f° r six 
simulations between low and high C a with 
M — 1 and 100 on a lattice of 40 x 60 and 
25 x 35 nodes, respectively. AT all (x) was cal- 
culated by including all valley sizes and the 
hierarchical smaller ones within a large val- 
ley (see Fig. 5). To obtain reliable average 
quantities we did 10 to 20 different simula- 
tions at each C a . In order to calculate the 
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Figure 6: The hierarchical valley size distri- 
butions A^ a n(x), for six simulations between 
low and high C a with M = 1 (O , □ , O) and 
M = 100 (A , < , V). The slope of the solid 
line is —1.9. 



valley sizes at large C a , we subtract the av- 
erage drift in the pressure signal due to vis- 
cous forces such that the pressure becomes a 
function that fluctuates around some mean 
pressure. 

By assuming a power law N a n(x) oc x~ Ta11 
our best estimate from Fig. 6 is r all = 1.9 ± 
0.1, indicated by the slope of the solid line. 
At low x m Fig- 6, typically only one tube 
is invaded during the valley and we do not 
expect the power law to be valid. Similar 
results were obtained when calculating the 
hierarchical distribution of the time periods 
T of the valleys, denoted as N a n(T). 

In invasion percolation (IP) the distribu- 
tion of burst sizes N(s), where s denotes 
the burst size, is found to obey the scaling 
relation [27,28,39,40] 

N(s)KS-'g(s°(f -f c )). (11) 



Here f c is the percolation threshold of the 
system and g(x) is some scaling function, 
which decays exponentially when x 3> 1 
and is a constant when x — > 0. r' is re- 
lated to percolation exponents like r' = 
1 + Df/D - l/(Du) [40], where D f and D 
is the fractal dimension of the front and the 
mass of the percolation cluster, respectively. 
v is the correlation length exponent in per- 
colation theory and a = l/{vD) [41]. In 
Eq. (11) a burst is defined as the connected 
structure of sites that is invaded following 
one root site of random number /o, along 
the invasion front. All sites in the burst 
have random numbers smaller than / , and 
the burst stops when the random number 
of the next site to be invaded is larger than 
fo [42]. 

By integrating Eq. (11) over all fo in the 
interval [0, f c ] Maslov [43] deduced a scal- 
ing relation for the hierarchical burst size 
distribution iV a ii(s) following 

/V a ii(s) oc (12) 

where r a u = 2. 

In the low C a regime in Fig. 6, the dis- 
placements are in the capillary dominated 
regime and the invading fluid generates a 
growing cluster similar to IP [4,7-9]. In 
this regime we also have that x s [28] 
and hence N^x) corresponds to iV a n(s) in 
Eq. (12). Thus, in the low C a regime we ex- 
pect that iVaji(x) follows a power law with 
exponent r a u = 2 which is confirmed by our 
numerical results. Similar results were ob- 
tained in Ref. [28]. 

The evidence in Fig. 6, that r a n does not 
seem to depend on C a , is very interesting. 
At high C a when M = 0.01 an unstable vis- 
cous fingering structure generates and when 
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M > 1 a stable front develops. It is an open 
question how these displacement processes 
map to the proposed scaling in Eq. (12). We 
note that in the high C a regime the relation 
X oc s may not be correct and T is preferred 
when computing 7V a n. However, the simu- 
lations show that Ngi\(x) ~ iV a ii(T) even at 
high C a . 

In Ref. [43] r a n was pointed out to be 
super universal for a broad class of self- 
organized critical models including IP. The 
result in Fig. 6 indicates that the simulated 
displacements might belong to the same su- 
per universality class even at high injection 
rates where there is no clear mapping be- 
tween the displacement process and IP. 

Basaket al. [44,45] performed four drain- 
age experiments where they used a 110 x 180 
mm transparent porous model consisting 
of a mono-layer of randomly placed glass 
beads of 1 mm, sandwiched between two 
Plexiglas plates. The experimental setup 
was similar to the one used in Ref. [27]. 
The model was initially filled with a water- 
glycerol mixture of viscosity 0.17 P. The 
water-glycerol mixture was withdrawn from 
one of the short side of the system at con- 
stant rate by letting air enter the system 
from the other short side. The pressure 
in the water-glycerol mixture on the with- 
drawn side was measured with a pressure 
sensor of our own construction. 

From the recorded pressure signal we cal- 
culated the hierarchical distribution of time 
periods of the valleys, iV all (T). At low C a 
this corresponds to iV all (s) in Eq. (12). Be- 
cause of the relative long response time of 
the pressure sensor, rapid and small pres- 
sure jumps due to small bursts are presum- 
ably smeared out by the sensor and the 
recorded pressure jumps are only reliable 
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Figure 7: The hierarchical distribution 
A^n (T) of the valley time T during a burst 
for experiments (open symbols) and simu- 
lations (filled symbols) at various C a with 
M = 0.017 and M = 0.01, respectively. 
The slope of the solid line is —1.9. 



for larger bursts. Hence, from the recorded 
pressure signal T appears to be a better es- 
timate of the burst sizes than \- 

In Fig. 7 we have plotted the logarithm 
of iV a ii(T) for experiments (open symbols) 
and simulations (filled symbols) performed 
at four different C a , respectively. To col- 
lapse the data N^T) and T were nor- 
malized by their means. In the simula- 
tions M = 0.01 while in the experiments 
M = 0.017 where we have assumed air 
to have viscosity 0.29 xl0~ 2 P. We observe 
that the experimental result is consistent 
with our simulations and we conclude that 
N all (T) oc T 19±01 . This confirms the scal- 
ing of N a ii(x) in Fig. 6. 

From Figs. 6 and 7 we conclude that 
r a n = 1.9 ± 0.1 for all displacement sim- 
ulations going from low to high injection 
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rates when M = 0.01, 1, and 100. This 
is also confirmed by drainage experiments 
performed at various injection rates with 
M = 0.017. The evidence that r a n is in- 
dependent of the injection rate, may indi- 
cate that the displacement process belongs 
to the same super universality class as the 
self-organized critical models in Ref. [43] 
= 2), even where there is no clear map- 
ping between the displacement process and 
IP. 

4 STABILIZATION 

MECHANISMS OF 
THE FRONT 

When the displacements are oriented out of 
the horizontal plane and the injection rate is 
low, gravity acting on the system, may sta- 
bilize the front due to density differences be- 
tween the fluids. Several authors [30,46-48] 
have confirmed, by experiments and sim- 
ulations, that the saturated front width 
w s scales with the strength of gravity like 
w s oc B ~ v ^ 1+V \ Here B Q (Bond number) 
is the ratio between gravitational and capil- 
lary forces, given by B Q = Apga 2 /^, where 
Ap is the density difference between the flu- 
ids, g the acceleration due to gravity, a the 
average pore size, and 7 the fluid-fluid in- 
terface tension. Furthermore, v denotes the 
correlation length exponent in percolation. 

A similar consensus concerning the sta- 
bilization mechanisms when the displace- 
ments are oriented within the horizontal 
plane, has not yet been reached. In hori- 
zontal displacements, viscous forces replace 
gravitational forces, and in the literature 
there exist different suggestions about the 



scaling of w s as function of C a . The capil- 
lary number C a is the ratio between viscous 
and capillary forces according to the defi- 
nition in Eq. (6). In 3D, where trapping 
of wetting fluid is assumed to be of little 
importance, Wilkinson [30] was the first to 
use percolation to deduce a power law like 
w s oc C a ~ a where a = vj{\ + t — (3 + v). 
Here t and f3 is the conductivity and order 
parameter exponent in percolation, respec- 
tively. Later, Blunt et al. [32] suggested in 
3D that a = v / (1 + 1 + v). This is identi- 
cal to the result of Lenormand [31] finding 
a power law as function of system size for 
the domain boundary in the C a -M plane 
between capillary fingering and stable dis- 
placement in 2D porous media. 

More recently, Xu et al. [33] used the ar- 
guments of Gouyet et al. [49] and Wilkin- 
son [30] to show that the pressure drop 
AP nw over a length Ah in the nonwetting 

phase of the front should scale as AP nw cx 
Ah t/v+D E -i-f3/v ( gee Fig gy Rere De de _ 

notes the Euclidean dimension of the space 
in which the front is embedded, i.e. in our 
case D E = 2. The pressure drop in the wet- 
ting phase AP W , was argued to be linearly 
dependent on Ah due to the compact phase 
there (Fig. 8). In Ref. [32] Blunt et al. also 
suggested a scaling relation for AP nw , how- 
ever, in 3D they found AP nw oc Ah t/u+1 . 
This is different from the result of Xu et al. 
when De = 3. 

In the next section, Sec. 4.1, we present 
an alternative view on the displacement 
process from those initiated by Wilkin- 
son [30], but include the evidence that non- 
wetting fluid flows in separate strands (see 
Fig. 9). The alternative view leads to an- 
other scaling of AP nw than the one Xu et 
al. [33] would predict in 2D, and we show 
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Figure 8: A schematic picture of the front 
that travels across the system from the in- 
let to the outlet. AP nw and AP W denote the 
viscous pressure drop going from A to B in 
the nonwetting and wetting phase, respec- 
tively. A and B are separated a distance 
Ah along the short side of the system. 



that it may influence a in the scaling be- 
tween w s and C a . The new scaling of AP nw 
is supported by numerical experiments us- 
ing the network model. 

4.1 LOOPLESS STRANDS 

Figure 9 shows two typical displacement 
structures that were obtained from simu- 
lations at low and high C a on a lattice of 
40 x 60 nodes with M — 1. From the fig- 
ure we observe that the nonwetting fluid 
(dark grey and black) generates patterns 
containing no closed loops. That means, 
following a path of nonwetting fluid will 
never bring us back to the starting point. 
The nonwetting fluid also flows in separate 
loopless strands, indicated as black tubes in 
Fig. 9. The loopless structures in Fig. 9 are 
a direct consequence of the evidence that 




Figure 9: Two displacement structures of 
simulations at high C a = 3.9 xlO -4 (top) 
and low C a = 1.6 xlO -5 (bottom) before 
breakthrough of nonwetting fluid. The lat- 
tice size is 40 x 60 nodes and M — 1. 
The nonwetting fluid (dark grey and black) 
is injected from below and wetting fluid 
(light grey) flows out along the top row. 
The black tubes denote the loopless strands 
where nonwetting fluid flows and the dark 
grey tubes indicate nonwetting fluid unable 
to flow (i.e. dead ends) due to trapped re- 
gions of wetting fluid. Note the few fluid 
supplying strands from the inlet to the 
frontal region at low C a compared to the 
case at high C a . 
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a tube filled with wetting fluid and sur- 
rounded on both sides by nonwetting fluid is 
trapped due to volume conservation of wet- 
ting fluid [50]. We note that this evidence 
may easily be generalized to 3D, and there- 
fore our arguments should apply there too. 
We also note that trapping of wetting fluid 
is more difficult in real porous media due 
to a more complex topology of pores and 
throats there. Loopless IP patterns have 
earlier been observed in Refs. [51-53]. 

From Fig. 9 we may separate the dis- 
placement patterns into two parts. One 
consisting of the frontal region continuously 
covering new tubes, and the other consist- 
ing of the more static structure behind the 
front. The frontal region is supplied by 
nonwetting fluid through a set of strands 
that connect the frontal region to the inlet. 
When the strands approach the frontal re- 
gion they are more likely to split. Since we 
are dealing with a square lattice, a splitting 
strand may create either two or three new 
strands. As the strands proceed upwards in 
Fig. 9, they split repeatedly until the frontal 
region is completely covered by nonwetting 
strands. 

On IP patterns without loops [51,53,54] 
the length / of the minimum path between 
two points separated an Euclidean distance 
R scales like I oc R Da where D s is the frac- 
tal dimension of the shortest path. We as- 
sume that the displacement pattern of the 
frontal region for length less than the cor- 
relation length (in our case w s ) is statis- 
tically equal to IP patterns in Ref. [51]. 
Therefore, the length of a strand in the 
frontal region is proportional to Ah Da when 
Ah is less than w s . If we assume that on 
the average every tube in the lattice has 
same mobility (k^/ /i^), this causes the fluid 



pressure within a single strand to drop like 
Ah Ds as long as the strand does not split. 
When the strand splits volume conserva- 
tion causes the volume fluxes through the 
new strands to be less than the flux in the 
strand before it splits. Hence, following a 
path where strands split will cause the pres- 
sure to drop as Ah K where k < D s . From 
Fig. 9, we note that at high C a the lengths 
of individual strands in the frontal region 
approach the minimum length due to the 
tubes. Therefore, in this limit finite size ef- 
fects are expected to cause D s — > 1. 

From the above arguments we conclude 
that the pressure drop AP nw , in the nonwet- 
ting phase of the frontal region (that is the 
strands) should scale as AP nw oc Ah K where 
k < D s . In 2D two different values for 
D s have been reported: D s = 1.22 [53, 54] 
for loopless IP with and without trapping 
and growing around a central seed, and 
D s = 1.14 [51] for the single strand con- 
necting the inlet to the outlet when non- 
wetting fluid percolates the system. In 
2D the result of Xu et al. [33] would give 
k = tjv + D E — 1 — P/u 1.9 where we 
have inserted t = 1.3, v = 4/3, (3 = 5/36, 
and .De = 2. This result for k is larger 
than D s and therefore not compatible with 
our arguments. 

To confirm numerically that k < D s , 
we have calculated the difference in cap- 
illary pressure AP C between menisci along 
the front by using our network model. AP C 
as function of Ah is calculated by taking 
the mean of the capillary pressure differ- 
ences between all pairs of menisci along the 
front, separated a distance Ah in the direc- 
tion of the displacement. AP C of a pair of 
menisci is calculated by taking the capillary 
pressure of the meniscus closest to the inlet 
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Figure 10: AP C as function of Ah for three 
different C a 's with M = 100 and 1 on lat- 
tices of 25 x 35 and 40 x 60 nodes, respec- 
tively. 



minus the capillary pressure of the meniscus 
closest to the outlet. 

Figure 10 shows AP C as function of Ah 
for simulations performed at high, interme- 
diate and low C a , with M = 1 or 100. The 
simulations with M = 100 were performed 
on a 25 x 35 nodes lattice with fx nw = 10 
P, /jl w = 0.10 P, and 7 = 30 dyn/cm. 
The disorder was introduced by choosing 
the tube radii at random in the interval 
0.05<i < Tjj < d. The tube length was 
d — 0.1 cm. The simulations with M — 1 
were performed on a distorted lattice of 
40 x 60 nodes where 0.02 cm < d - < 0.18 
cm and = dij/2a with a = 1.25. Here 
l^nw — — 0.5 P. We did 10-30 simu- 
lations at each C a to obtain reliable aver- 
age quantities. From the plot we observe 
that AP C increases roughly linearly as func- 
tion of Ah. The simulations also show that 
AP C ~ AP nw — AP w and that AP W depends 



linearly on Ah due to the compact wetting 
phase there (see Fig. 8). Hence our sim- 
ulations give that AP nw ~ AP C oc Ah K 
where k ~ 1.0. This supports the argu- 
ments giving k < D s , and we conclude that 
earlier proposed theories [30-33] which do 
not consider the evidence that nonwetting 
fluid flows in strands, are incompatible with 
drainage when strands are important. 

To save computation time and thereby be 
able to study AP C on larger lattices in the 
low C a regime, we have generated bond in- 
vasion percolation (IP) patterns with trap- 
ping on lattices of 200 x 300 nodes. The 
bonds in the IP lattice correspond to the 
tubes in the network model. The occupa- 
tion of bonds started at the bottom row, 
and the next bond to be occupied was 
always the bond with the lowest random 
number / from the set of empty bonds along 
the invasion front. We applied a small gra- 
dient in the random numbers of the bonds 
to stabilize the front [30,47]. 

When the IP patterns became well de- 
veloped with trapped (wetting) clusters of 
sizes between the bond length and the front 
width, the tubes in our network model were 
filled with nonwetting and wetting fluid ac- 
cording to occupied and empty bonds in 
the IP lattice. Furthermore, the radii r 
of the tubes were mapped properly to the 
random numbers / of the bonds by first re- 
moving the gradient in / and then assigning 
r — 1 — f [36]. The last transformation is 
necessary because in the IP algorithm the 
next bond to be invaded is the one with the 
lowest random number, opposite to the net- 
work model, where the widest tube will be 
invaded first. 

Having initialized the tube network, the 
network model was started and the simu- 
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lations were run a limited number of time 
steps while recording AP C . The number of 
time steps where chosen sufficiently large 
to let the menisci along the front adjust 
according to the viscous pressure set up 
by the injection rate. By this method we 
save the computation time that would have 
been required if the displacement patterns 
should have been generated by the network 
model instead of the much faster IP algo- 
rithm. However, to make this method self- 
consistent we have to assume that the IP 
patterns are statistically equal to the corre- 
sponding structures that would have been 
generated by the network model. 

Totally we generated four IP patterns 
with different sets of random numbers and 
every pattern was loaded into the network 
model. The result of the calculated AP C 
versus Ah is shown in Fig. 11 for C a = 
9.5 xl(T 5 and M = 100. If we assume 
a power law AP C oc Ah K , we find k = 
1.0 ± 0.1. The slope of the straight line in 
Fig. 11 is 1.0. We have also calculated AP C 
for C a = 2 x 10~ 6 with M = 1 and M = 100 
by using one of the generated IP patterns. 
The result of those simulations is consistent 
with Fig. 11 which indeed show the similar 
behavior as the results in Fig. 10. 

An important issue, arising at low C a , is 
the effect of bursts on the capillary pres- 
sure. A burst occurs when a meniscus along 
the front becomes unstable and nonwetting 
fluid abruptly covers new tubes [28]. The 
tube where the burst initiates will during 
the burst, experiences a much higher fluid 
transport relative to tubes far away. De- 
scribing the pressure behavior between the 
tube of the burst and the rest of the front is 
nontrivial. However, simulations show that 
even during bursts, we find that AP C in- 
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Figure 11: The logarithm of AP C as func- 
tion of the logarithm of Ah for simulations 
initiated on IP patterns at C a = 9.5 xl0~ 5 
and M = 100. The slope of the solid line is 
1.0. 



creases linearly with Ah. 

The evidence that k ~ 1.0 may influence 
the exponent a in «j s oc C a ~ a . Assum- 
ing Darcy flow where the pressure drop de- 
pends linearly on the injection rate, we con- 
jecture that AP C oc C a Ah K . Here AP C de- 
notes the capillary pressure difference over 
a height Ah when the front is stationary. 
That means, AP C excludes situations where 
nonwetting fluid rapidly invades new tubes 
due to local instabilities (i.e. bursts). The 
above conjecture is supported by simula- 
tions showing that in the low C a regime 
AP C oc C a Ah K where k ~ 1.0. Note, that 
AP C AP C in Fig. 10 since the latter in- 
cludes both stable situations and bursts. 

At sufficiently low C a where only the 
strength of the capillary pressure decides 
which tube should be invaded or not, we 
may map the displacement process to per- 
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eolation giving AP C oc f-f c oc f" 1 /" [30,47, 
49]. Here / is the random numbers in the 
percolation lattice, f c is the critical percola- 
tion threshold, and £ oc w s is the correlation 
length. Combining the above relations for 
AP C gives w s oc C a ~ a where a = i//(1+vk). 
In 2D v = 4/3 and by inserting k = 1.0 we 
obtain a ~ 0.57. Note that this is differ- 
ent to results suggested in Refs. [30,32,33] 
giving a « 0.37-0.38 in 2D. 

At high C a the nonwetting fluid is found 
to invade simultaneously everywhere along 
the front, and consequently the front never 
reaches a stationary state [36]. In this limit 
simulations show a nonlinear dependency 
between AP C and C a . Therefore, in the high 
C a regime it is not clear if the above map- 
ping to percolation is valid, and we expect 
another type of relation between w s and C a . 

4.2 COMPARISON WITH 
EXPERIMENTS 

Frette et al. [55] have performed two- 
phase drainage experiments in a 2D porous 
medium with viscosity matched fluids 
(M = 1). They reported on the stabi- 
lization of the front and measured the sat- 
urated front width w s , as function of C a . 
Their best estimate on the exponent in 
w s oc C a ~ a was a = 0.6 ± 0.2. This is 
consistent with the above conjecture (a = 
0.57), however, corresponding simulations 
on 40 x 60 nodes lattices give a = 0.3 ±0.1. 
Figure 12 contains the data of Frette et 
al. (filled circles) and the result of our 
simulations (open diamonds). The simu- 
lations are performed at C a > 1.0 xl0~ 5 
while most of the experiments where done 
at C a < l.OxlO" 5 . Since the range of the 
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Figure 12: log 10 (w s ) as function of log 10 (C a ) 
for experiments from [55] and simulations 
on the lattice of 40 x 60 nodes. In both 
experiments and simulations M — 1. The 
slope of the solid and dashed line is -0.6 and 
-0.3, respectively. 



two does not overlap it is difficult to com- 
pare the result of the simulations with those 
of the experiments. However, the change 
in a from 0.6 to 0.3, might be consistent 
with a crossover to another behavior at high 
C a according to the discussion in Sec. 4.1. 
We also note, that in the simulations at 
C a — 1.0 xlO -5 , the front width approaches 
the maximum width due to the system size, 
making it difficult to observe any possible 
a ~ 0.57 regime at low C a . We emphasize 
that more simulations on larger systems and 
at lower C a are needed before any conclu- 
sion on a can be drawn. 

4.3 DISCUSSION 

The evidence that the nonwetting fluid dis- 
places the wetting fluid in a set of loopless 
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strands opens new questions about the dis- 
placement process. Returning to Fig. 9 it 
is striking to observe the different patterns 
of strands at high and low C a . At low C a 
few strands are supplying the frontal region 
with nonwetting fluid, and the strands split 
many times before the whole front is cov- 
ered. At high C a the horizontal distance 
between each strand in the static structure 
is much shorter, and only a few splits are 
required to cover the front. We conjecture 
that the average horizontal distance be- 
tween the fluid supplying strands depends 
on the front width. However, further inves- 
tigation of the displacement patterns is re- 
quired before any conclusions can be drawn. 

So far the arguments in Sec. 4.1 only con- 
sider displacements where the nonwetting 
strands contain no loops. A very interesting 
question that has to be answered is: What 
happens to k when different strands in the 
front connect to generate loops. In ordinary 
bond or site percolation loops generally oc- 
cur. Loops are also observed in experiments 
corresponding to those of Frette et al. [55]. 
In the experiments it is more difficult to 
trap wetting fluid due to the more complex 
topology of pores and throats (see Fig. 1). 
Consequently, loops will more easily gen- 
erate there, than in the case of a regular 
square lattice. Loops might also be created 
when neighboring menisci along the front 
overlap and coalesce depending on the wet- 
ting properties of the nonwetting fluid [5,6]. 

As a first approximation we conjecture 
that creation of loops will not cause k to 
change significantly. Note that in the front 
the different nonwetting strands connecting 
to each other to create loops, must at some 
later time split. Otherwise successive con- 
nections will cause the different strands to 



coalesce into one single strand of nonwet- 
ting fluid. Moreover, after the front width 
has saturated, the number of places where 
different strands connect must on average 
be balanced by the number of places where 
strands split. Therefore, we believe that 
the influence on k due to connections (i.e. 
loops) will be compensated by the splits and 
the overall behavior of k will remain the 
same. We emphasize that further simula- 
tions and experiments are required to in- 
vestigate the effect of loops on k. 

According to the discussion in Sec. 4.3, 
the evidence that the displacement patterns 
consist of loopless strands may easily be 
generalized to 3D. Therefore our arguments 
giving k < D s , might be valid in 3D as 
well. Note also that in 3D it is less proba- 
ble that different strands meet. Hence, even 
if they were supposed to connect to create 
loops, the number of created loops are ex- 
pected to be few. In 3D the fractal dimen- 
sion of the shortest path for loopless IP is 
D s = 1.42 [53]. 

5 CONCLUSION 

We conclude that our 2D network model 
properly simulates the temporal evolution 
of the pressure in the fluids during drainage. 
We have seen that capillary forces situ- 
ated around isolated and trapped regions of 
wetting fluid, contribute to the total pres- 
sure across the lattice, as well as the capil- 
lary fluctuations due invasion of nonwetting 
fluid along the front. 

We have found that the model reproduces 
the typical burst dynamics at low injection 
rates and we have investigated the statis- 
tics of the bursts by calculating the hierar- 
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chical valley size distribution N a n(x)- We 
conclude that Ng^(x) follows a power law, 
^aii(x) oc x~™ where r aU = 1.9 ± 0.1 is 
independent of the injection rate and vis- 
cosity ratio. Similar result is obtained from 
experiments. At low injection rates the 
result is consistent with the prediction in 
Ref. [43] (r a n = 2), which was deduced for 
a broad spectrum of different self-organized 
critical models including IP. The evidence 
that r a u is independent of the injection rate, 
may indicate that the displacement pro- 
cess belongs to the same super universality 
class as the self-organized critical models in 
Ref. [43], even where there is no clear map- 
ping between the displacement process and 
IP. 

We have also simulated the behavior of 
the capillary pressure along the front. Sim- 
ulations show that the capillary pressure 
difference AP C between two points along 
the front varies almost linearly as function 
of distance Ah in the direction of the dis- 
placement. The numerical result supports 
arguments based on the observation that 
nonwetting fluid flows in separate strands 
where wetting fluid is displaced. From the 
arguments we find that AP C oc Ah K where 
k < D s . Here D s denotes the fractal dimen- 
sion of the nonwetting strands. 

Several attempts have been made to 
describe the stabilization mechanisms in 
drainage due to viscous forces, however, 
none of them consider the evidence that 
nonwetting fluid displaces wetting fluid 
through strands. Therefore, we conclude 
that earlier suggested theories fail to de- 
scribe the stabilization of the invasion front 
when strands dominate the displacements. 
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